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Abstract We investigate the dynamic mechanisms of generation of subthreshold and 
phase resonance in two-dimensional linear and linearized biophysical (conductance- 
based) models, and we extend our analysis to account for the effect of simple, but not 
necessarily weak, types of nonlinearities. Subthreshold resonance refers to the ability 
of neurons to exhibit a peak in their voltage amplitude response to oscillatory input 
currents at a preferred non-zero (resonant) frequency. Phase-resonance refers to the 
ability of neurons to exhibit a zero-phase (or zero-phase- shift) response to oscillatory 
input currents at a non-zero (phase-resonant) frequency. We adapt the classical phase- 
plane analysis approach to account for the dynamic effects of oscillatory inputs and 
develop a tool, the envelope-plane diagrams, that captures the role that conductances 
and time scales play in amplifying the voltage response at the resonant frequency 
band as compared to smaller and larger frequencies. We use envelope-plane diagrams 
in our analysis. We explain why the resonance phenomena do not necessarily arise 
from the presence of imaginary eigenvalues at rest, but rather they emerge from the 
interplay of the intrinsic and input time scales. We further explain why an increase 
in the time-scale separation causes an amplification of the voltage response in addi- 
tion to shifting the resonant and phase-resonant frequencies. This is of fundamental 
importance for neural models since neurons typically exhibit a strong separation of 
time scales. We extend this approach to explain the effects of nonlinearities on both 
resonance and phase-resonance. We demonstrate that nonlinearities in the voltage 
equation cause amplifications of the voltage response and shifts in the resonant and 
phase-resonant frequencies that are not predicted by the corresponding linearized 
model. The differences between the nonlinear response and the linear prediction in- 
crease with increasing levels of the time scale separation between the voltage and the 
gating variable, and they almost disappear when both equations evolve at comparable 
rates. In contrast, voltage responses are almost insensitive to nonlinearities located in 
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the gating variable equation. The method we develop provides a framework for the 
investigation of the preferred frequency responses in three-dimensional and nonlinear 
neuronal models as well as simple models of coupled neurons. 

1 Introduction 

Rhythmic oscillations have been observed in various areas of the brain and have been 
implicated in cognition and motor behavior [1-4] in both health and disease [5]. 
Network oscillations result from the cooperative activity of the participating neurons 
[3]. Many neuron types possess membrane potential oscillatory properties [4], which 
emerge either in the form of intrinsic subthreshold oscillations (STOs) [4, 6-9], or 
subthreshold resonance [10-20], or both [12, 15, 21, 22]. Subthreshold resonance 
refers to the ability of neurons to exhibit a peak in their voltage amplitude response to 
oscillatory input currents at a preferred non-zero (resonant) frequency [10]. Intrinsic 
STOs in isolated neurons emerge spontaneously or in response to tonic (DC) current 
inputs, and primarily reflect interactions among the neuron's intrinsic currents. In 
contrast, subthreshold resonance results from the interaction between these intrinsic 
currents and an oscillatory input currents. Because of this, subthreshold resonance 
has been implicated in the generation of oscillations at the network level [23-27]. 

The relation between intrinsic STOs and subthreshold resonance is still an open 
question. For some neuron types, STOs and subthreshold resonance have been shown 
to result from the same mechanism [12]. However, theoretical and experimental stud- 
ies have demonstrated that they are not equivalent phenomena [1 1, 28]: neurons may 
exhibit one and not the other [10, 11, 21]. Furthermore, standard calculations for 
linear models show that their natural (intrinsic) and resonant frequencies do not gen- 
erally coincide except in some, rather restricted, parameter regimes [1 1, 29] (see also 
our discussion in Appendix A.3). 

The phase-shift (or phase) of the neuronal voltage response to subthreshold os- 
cillatory input currents has received less attention that the corresponding amplitude 
response [11, 14, 29]. This despite the fact that phases are expected to play a ma- 
jor role in determining the synchronization properties of neuronal networks [30]. 
A zero-phase response indicates that both voltage output and current input peak 
at the same time, thus generating in-phase synchronized patterns. We use the term 
phase-resonance to refer to the ability of neurons to exhibit a zero-phase response to 
oscillatory inputs at a non-zero (phase-resonant) frequency. The resonant and phase- 
resonant frequencies do not generally coincide [29] for neuronal models (they do so 
for circuits in parallel but not for circuits in series as neuronal models are). In ad- 
dition, resonance may occur in the absence of phase-resonance [29] (see also our 
discussion in Appendix A.3). 

The properties of subthreshold resonance have been investigated in many sys- 
tems [9, 11, 14, 15, 20, 21, 25, 29]. Theoretical studies have focused on simulations 
of conductance-based models and the analysis of the impedance profiles (curves of 
amplitude and phase- shift as a function of the input frequency) for the correspond- 
ing linearized models around resting potential [11, 14-16, 18, 31-35]. However, the 
mechanisms underlying the generation of resonance and phase-resonance in neurons 
are not fully understood. This is partly because adequate tools are lacking. 
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These mechanisms can be addressed from two different, but complementary per- 
spectives: biophysical and dynamic. The former focuses on the role of the ionic cur- 
rents and their biophysical properties in shaping the neuron's voltage response, and 
it has been discussed in terms of the so-called resonant and amplifying ionic cur- 
rents (see Sect. 2.5) [10, 11]. The latter focuses on (i) the geometric properties of the 
neuronal models in terms of the nullclines and phase planes, and (ii) the interaction 
between the neuron's intrinsic time scales and the time scales associated with the 
input currents to produce optimal voltage responses in both amplitude and phase. 

In [29] we have identified the basic biophysical mechanisms of generation of res- 
onance and phase-resonance in two-dimensional linear and linearized conductance- 
based models, and we have conducted a thorough study of the properties of the volt- 
age response in terms of the biophysical parameters. In particular, we have shown 
how changes in the maximal conductances affect the resonant and phase-resonant 
frequencies and other attributes of the impedance amplitude and phase profiles in 
ways that are not always intuitive. 

The goal of this paper is to investigate the dynamic mechanisms that give rise to 
resonance and phase-resonance. Specifically, we aim to identify the basic dynamic 
and geometric principles that govern the generation of these phenomena in order to 
understand (i) how the resonant and phase-resonant frequencies are selected, (ii) how 
the voltage response is amplified at the resonant frequency band, (iii) how these prop- 
erties are affected by changes in parameters (e.g., maximal conductances, time con- 
stants), and (iv) how resonance and phase-resonance are related to intrinsic STOs. 

We use dynamical systems tools and numerical simulations, and we extend the 
classical phase-plane analysis approach to generate a tool (envelope-plane diagrams) 
that enables us to address the mechanistic issues described above. For simplicity, in 
this paper we focus on two-dimensional linear systems and we illustrate how the ideas 
we develop here can be used to investigate linearized conductance-based models and 
nonlinear models. 

From a dynamic perspective, we view the properties of the forced oscillations as 
reflecting the interaction between the forced neuron's intrinsic time scales, which 
result from a combination of its biophysical properties and the time scales of the 
input. The time scales that emerge from these interactions determine the amplitude 
and phase of the voltage response for each value of the input frequency. The resonant 
and phase-resonant frequencies correspond to the emergent time scales that allow the 
neuron to maximize its voltage amplitude response and to peak in phase with the 
input, respectively. Therefore, we aim to elucidate how these emergent time scales 
are generated, how intrinsic and emergent time scales are related, and how they are 
affected by changes in the model parameter. 

It is clear that for linear systems the voltage response can be computed ana- 
lytically. However, due to its complexity, the dependence of the voltage response 
properties with the model parameters is not straightforward, not even for linear two- 
dimensional systems [29] . Additionally, due to the same complexity, it is difficult to 
extract from the closed-form solutions a mechanistic understanding of the neuron's 
resonant properties in terms of the time scales. The use of geometric tools aids in this 
effort. Since the phase plane contains information about the biophysical properties of 
neurons through the structure of the nullclines and time scales, dynamical systems 
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tools allow us to make statements that are valid for generic classes of biophysical 
models. 



2 Methods 

2.1 Conductance-Based Models 

We consider conductance-based models of Hodgkin-Huxley type [36]. The current- 
balance equation is given by 

= -/ L -h-h + 4pp + /in(0, (1) 

where V is the membrane potential (mV), t is time measured in msec, C is the mem- 
brane capacitance (uF/cm 2 ), 7 app is the applied bias (DC) current (uA/cm 2 ), I m (t) is 
a time-dependent input current (uA/cm 2 ), 7l = G^(V — £"l) is the leak current, and 
Ij 0 = 1, 2) are ionic currents of the form 

Ij = GjXj(V-Ej) (2) 

with maximal conductance Gj (mS/cm 2 ) and reversal potentials Ej (mV), respec- 
tively. The ionic currents (2) we consider here are restricted to have a single gating 
variable xj and to be linear in xj . This choice is motivated by the persistent sodium 
(^NapX h- (7h) and M-type (7m) currents found in several neurons that exhibit sub- 
threshold resonance [15, 37-39] (see Appendix B). All gating variables x obey a first 
order differential equation of the form 

dx_ _ Xqo(V)-x 
dt ~ r x (V) ' 



where x^iV) and r x (V) are the voltage-dependent activation/inactivation curves and 
time scales, respectively. For external sinusoidal inputs we use the following notation: 

27tf 

knit) = A in sin(tf 0 with Q = -y-, (4) 

where T = 1000 msec and [/] = Hz. 

In this paper we focus on two-dimensional models having one dynamic gating 
variable (x\) and, possibly, an additional gating variable evolving on a fast time scale 
for which the adiabatic approximation X2 = ^2,00 (V) is made. Additional fast cur- 
rents can be including without significantly changing the formalism. The investiga- 
tion of three-dimensional systems is beyond the scope of this paper. 

2.2 Linearized Conductance-Based Models 

We follow Richardson et al. [11] and linearize the autonomous part of system (l)-(3) 
(/ in = 0) around the fixed point (V , x\) by defining 

v = V — V and w = — — (5) 
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where x\ = x\, OQ (V). The linearized equations are [11] 

dv 

C— = -g L V ~ g\W + h n (t), (6) 

dt 

- dW ,H\ 

t 1 — = v-w, (7) 
dt 

where the effective conductances #l and g\, and time constant f\ are defined by 

g L = G L + Gix hoo (V) + G 2 x 2 ,oo(V) + g 2 , (8) 

gj=G j (V-Ej)x' Jt „<y), j = 1,2 (9) 

and 

fi = n(V). (10) 

The sign of the effective ionic conductances gj ( j = 1 , 2) determines whether 
their associated gating variables are resonant (gj > 0) or amplifying (gj < 0) [10, 
11]. The effective conductance #l includes information not only about the original 
leak conductance Gl but also about the voltage-dependent ("passive") terms of the 
ionic currents. The ionic currents I\ and h each contribute a positive term to The 
ionic current h contributes to #l with an additional term (#2) due to its instantaneous 
dynamics. If x 2 is amplifying and g2 < 0 is large enough in absolute value, then #l 
is negative. Note that the gating variable w in (5) has units of voltage. 

2.3 Rescaled Linearized Models 

We rescale system (6)-(7) to aid in the analysis and to reduce the number of parame- 
ters that effectively govern its dynamics (without loss of information). The rescaling 
we use here focuses on the geometric properties of the system and the time- scale 
separation between the participating variables, and is amenable to analysis using dy- 
namical systems tools (phase-plane analysis). A different scaling, appropriate for ad- 
dressing biophysically related questions, has been used in [1 1, 29]. 

We define the following dimensionless time and (dimensional) voltage variables: 

t = — t, v = — v = — , (11) 
C gi a 



parameters 



and 



81 C nn\ 
a = — , € = - , (12) 



I m (t) = A in sin(2n fi/T) with 

7 ~ Tg L , 2nf ( 13 > 

Ai n = , T = and Q = . 

gl C T 
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a b 




G h G 
h q 

Fig. 1 Dimensionless parameters a and 6 as a function of the biophysical resonant conductances in two 
representative conductance-based models. The models are described in Appendix B. a 1^ + /Nap model 
with slow /h and fast /Nap (biophysical conductances: G^ and G p , respectively), b /k s + ^Nap model with 
slow /ks and fast /Nap (biophysical conductances: G q and G p , respectively). For small enough values of 
Gp (dashed curves) both a and e are positive for all values of Gh and G q (a and b, respectively), while 
for larger values of G p (solid curves) both a and e are negative for large values of G^ and small values 
of Gq (in a and b, respectively). Note that although /^ and /k s are resonant currents, both a and e exhibit 
different monotonic properties as Gh (a) and G q (b) increase 



Substituting into (6)-(7) and dropping the "hat" signs we obtain 

dv 

— = -v - w + A in sm(Qt), (14) 
at 



— = e\av — w]. (15) 
dt 

Geometrically, the parameter a is the slope of the w-nullcline and can be thought 
of as representing the strength of the gain of the feedback in the linearized system. 
The parameter € represents the time-scale separation between v and w. 

Since for resonant gating variables g\ > 0, the sign of both a and € depends on 
whether #l is positive or negative. In the absence of fast amplifying currents (G2 = 
g2 = 0), gL > 0 and then both a > 0 and € > 0. When an amplifying current is present 
and its contribution to #l is small enough, the sign of both a and € remains positive 
(Fig. 1, dashed curves). However, when stronger contributions of the fast amplifying 
current causes #l to be negative, the sign of both a and € are also negative (Fig. 1, 
solid curves). Resonance occurs in both cases [11, 29]. Since resonance becomes 
amplified as #l decreases [29], we expect resonance to be more amplified for negative 
values of both € and a as compared to positive ones. The cases including values of 
a and € having different signs are excluded from this study since the underlying 
autonomous system is either unstable (saddle) or stable (node) but does not exhibit 
resonance (Fig. 3). In both cases x\ is amplifying (g\ < 0). The case € < 0 and a > 0 
requires that both gh < 0 and gi < 0, while the case € > 0 and a < 0 requires that 
g L > 0 and gi < 0. 

2.4 Impedance and Impedance-Like Functions 

The voltage response of a linear system receiving sinusoidal current inputs of the 
form (4) is given by 

VoutC* ; /) = Aou t (/) sin(^ - 0CO), (16) 
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Fig. 2 Schematic diagrams of the impedance (a) and phase (b) profiles (impedance and phase as a 
function of the input frequency /). al Band-pass filter (resonance). a2 Low-pass filter (no resonance), 
bl Zero-frequency phase crossing (phase-resonance). b2 Monotonically increasing and positive phase (no 
phase-resonance), a The resonant frequency / res is the input frequency / at which the impedance Z(f) 
reaches its maximum Z max . The resonance amplitude Qz = Z max — Z(0) measures the resonance power. 
The half-width frequency band A\j2 is the length of the frequency interval in between / res and the input 
frequency value at which Z(/) = Z max /2, and measures the system's selectivity to incoming frequencies 
close to / re s • b The phase-resonant frequency /phas is the zero-crossing phase frequency. The minimum 
phase 0 m i n measures the magnitude of the negative phase 



where A out is the amplitude and 0 is the phase- shift (or phase), defined as the 
difference between the peaks of the current input I m (t\ f) and the voltage output 
Voat(f;f). 

Linear systems exhibit resonance if there is a peak in the amplitude of the 
impedance function Z(f) given by 

\ Z (f)\=^n. d7) 

^in 

at some positive (resonant) frequency / res . In what follows, we will refer to 
impedance amplitude simply as the impedance Z(f). In Fig. 2a we show repre- 
sentative graphs of the impedance function Z(f) for a model that does (panel al) 
and does not (panel a2) exhibit resonance. We characterize the impedance profiles 
using four parameters: (i) the resonant frequency / res , (ii) the maximum impedance 
Zmax = Z(/ res ), (iii) the resonance amplitude Qz = Z max — Z(0), and (iv) the half- 
width frequency band A 1/2, defined as the frequency interval in between / res and the 
input frequency value at which Z(f) = Z max /2. A 1/2 is a measure of the frequency 
selectivity. Neurons have a higher selectivity to inputs with frequencies around / res 
the sharper the graph of Z(f). In Fig. 2b we show two representative graphs of the 
phase 0 (/) where 0 vanishes at a non-zero value of / (panel bl) and 0 is always pos- 
itive (panel b2). We refer to the ability of the neuron to exhibit a zero-phase frequency 
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response at a non-zero frequency as phase-resonance and to the corresponding fre- 
quency as the phase-resonant frequency / p has- in panel bl, / p h as > 0- The voltage 
response is "advanced" and "delayed" for lower and higher frequency inputs, respec- 
tively. Although phase advance and phase delay are ambiguous concepts to describe 
phase differences between inputs and outputs in oscillatory systems, we still use them 
since typical phase differences lie in the range (— tt/2, tt/2). In panel b2, / p h as = 0, 
that is, the voltage response is delayed for all values of /. We characterize the phase 
profiles using two parameters: (i) the phase-resonant frequency / p has, and (ii) the 
minimum phase 0 m in- 

For nonlinear systems, or for linear systems with non- sinusoidal inputs, (17) does 
no longer provide an appropriate definition of the impedance function. Here we as- 
sume that the voltage output is periodic and has the same frequency as the input and 
we use the following definition: 

ry , n\ ^max(/) — ^min(/) /io\ 

L \J) = tt: > Uo) 

Z/ii n 

where V max (f) and Vmin(/) are the maximum and minimum of the oscillatory volt- 
age response V out (f) for each value of the input frequency /. For linear systems 
receiving sinusoidal inputs, (18) and (17) are equivalent. The resonant frequency / res 
is the peak frequency of Z(f) in (18). Similarly to the linear case the phase 

is computed as the distance between the peaks of the output and input normalized 
by the period. Note that (18) can be thought of as a filtered version of the impedance 
function computed using the so-called ZAP functions [10, 1 1, 28] that sweep through 
a given range of frequencies continuously over time. 



2.5 Resonant and Amplifying Ionic Currents 



Biophysically, subthreshold resonance has been argued to result from a combination 
of low- and high-pass filter mechanisms that have been described in terms of neural 
currents [10]. RC circuits act as low-pass filters. (As the input frequency increases the 
voltage amplitude response of passive neurons decreases from its resistance value, for 
zero input frequency, to zero, for input frequencies approaching infinity.) 

Ionic currents, or more precisely their associated gating variables, have been clas- 
sified into resonant (g\ > 0) and amplifying (g\ < 0) [10, 11]. Resonant gating vari- 
ables (e.g., hyperpolarization-activated h-currents [15-17] and slow potassium, M- 
type currents [14]) have the ability to create resonance by opposing voltage changes 
(negative feedback). Amplifying gating variables (e.g., persistent sodium currents 
[14-17] and high-threshold calcium currents [10]) generate a positive feedback ef- 
fect that enhances voltage changes but they do not create resonance [10]. Some cur- 
rents such as the low-threshold calcium current Ij have both resonant and amplifying 
gating variables [10]. 

The interaction between resonant and amplifying currents in shaping the neuronal 
voltage response (impedance and phase profiles) to oscillatory input currents is com- 
plex [29], sometimes non-intuitive, and involves not only changes in the maximal 
amplitude of the voltage response but also in other attributes including the resonant 
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and phase resonant frequencies. The role that different types of resonant and amplify- 
ing currents play in shaping the impedance profile has been recently clarified in [29] . 
An important outcome of this study is that the standard classification "resonant vs. 
amplifying" does not capture in its entirety the effect that changes in their biophysi- 
cal parameters have on the shapes of the impedance and phase profiles. For instance, 
while hyperpolarization-activated (h- or 7 n ) and M-type slow potassium (Iks) cur- 
rents have qualitatively similar effects on the impedance profile as the corresponding 
ionic conductances (G n and Gk s ), they may have opposite effects in the presence of 
a persistent sodium current (amplifying). Specifically, in the latter case, an increase 
in G n leads to an amplification of the voltage response, while an increase in Gks 
causes an attenuation of the voltage response. Additionally, an increase in the time 
constant associated to the gating variable, which increases the system's time-scale 
separation, does not only affect the resonant and phase-resonant frequencies, but it 
may also produce a significant amplification of the voltage response. 



3 Results 

3.1 Stability and Resonant Properties in the a-€ Rescaled System 

Here we review the stability and resonant properties of the rescaled equations (14)- 
(15) in terms of the parameters a and € for later use. These properties have been 
discussed in [11, 29] for a different rescaling based on dimensionless effective con- 
ductances. 

3.1.1 Stability Properties for the Autonomous System and Intrinsic Oscillations 

We first consider system (14)— (15) with A[ n = 0. The fixed point is given by (v, w) = 
(0, 0) and the eigenvalues are given by (Appendix A) 



-(1 + €) ± 7(1 -O 2 - 4c*6 

n,2 = 2 ' ^ ^ 

The stability diagram in the a-€ parameter space is presented in Fig. 3al. From (19), 
the fixed point is a focus if (1 — e) 2 — 4ae <0 and a node otherwise. Foci are stable 
if 1 + e > 0. Nodes are stable if 1 + € > 0 and e(l + a) > 0. If e(l + a) < 0 the fixed 
point is a saddle. The natural frequency of the damped oscillations (for stable foci) is 
given by 



T J4ae - (1 - 6) 2 

/nat=— tl With At = ^ £ -. (20) 

2tt 2 
3.1.2 Impedance, Phase, Resonance, and Phase-Resonance 

The impedance function (17) for system (14)— (15) with A- m > 0 is given by (Ap- 
pendix A) 

z m ~ [€(i+a)-^ 2 ] 2 + a + o 2 ^ 2 ' (21) 
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al a2 a3 




0 2 4 6 80 2 4 6 8 



Fig. 3 Stability and resonant properties for the reduced two-dimensional linear system (14)— (15). a Stabil- 
ity and resonance diagrams in the a-e parameter space, al Stability diagram for the autonomous system. 
The blue curves separate between regions with different stability properties. a2 Resonance diagram. The 
red curves separate between regions where the system does (above and below) and does not (middle) ex- 
hibit resonance. a3 Superimposed stability (blue curves) and resonance (red curves) diagrams showing 
that intrinsic oscillations and resonance may occur in the absence of the other, b Natural (/nat), resonant 
(/res), and phase-resonant (/ p has) frequencies as a function of a (bl) and e (b2) illustrating that these 
characteristic frequencies have different values (bl and b2) and different monotonic properties (b2) as the 
model parameters change, c Maximum impedance (Z max ) and resonance amplitude (Qz) as a function 
of a (cl) and e (c2) illustrating the two basic mechanisms of generation of resonance in 2D linear sys- 
tems, cl As a increases, resonance results from a decrease in both Z(0) and Z max with Z(0) decreasing 
faster than Z max . c2 As e decreases (time-scale separation increases), resonance results from an increase 
in Z m ax 

with Z(0) fixed 



The resonant frequency, if it exists, is given by 

/res = ^res with ^ res = J -€ 2 + y/ € 2 a(a + 2e + 2). (22) 
2tt v 

System (14)— (15) exhibits resonance if £2 res > 0. In order £? res to be defined, a(a + 
2£ + 2) > 0 and ^ € 2 a(a + 2^+2) > € 2 . The first condition is satisfied either if 

a > -2(1+0 (a>0) or a < -2(1 + 0 (a < 0). (23) 
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For the second condition to be satisfied, 

|a + l+€| >y2e 2 + 2e + l. (24) 

This yields either 

a > -1 - € + ^2e 2 + 2e + 1 or a < -1 - € -^2e 2 + 2€ + 1. (25) 

These regions are illustrated in the resonance diagram in Fig. 3a2. Figure 3a3 illus- 
trates that resonance may occur in the absence of intrinsic oscillations and vice versa. 
The phase 0 for system (14)— (15) is given by 

(j){Q) = tan" 1 ^ . (26) 

The phase-resonant frequency, if it exists, is given by 

T , 

/phas = — %ias with ^ phas = ^e(a-e). (27) 
Lit 

The natural, resonant, and phase-resonant frequencies do not necessarily coincide 
(Fig. 3b). In fact, they rarely do so. 

Figure 3c illustrates the two basic mechanisms of generation of resonance for 2D 
linear systems [29] in terms of a and € : (i) as a increases, resonance emerges from 
a combined and unbalanced decrease in both Z(0) and Z max with Z(0) decreasing 
faster than Z max (panel cl), and (ii) as € decreases (time-scale separation between 
v and w\ increases), resonance emerges from an increase in Z max with Z(0) un- 
changed. 



3.2 The Structure of the "Oscillatory Phase Plane" 



Here we consider system (14)— (15) with A- m > 0. For the type of analysis we present 
in this paper it is useful to rescale time by defining 

i=ft (28) 

in order to separate the effect of the input frequency / from the input's time depen- 
dence. Substituting (28) into (14)— (15), and dropping the "hat" sign from the time t, 
we obtain 

1 

-w + A t ], (29) 
- w], (30) 

ai j 
where the sinusoidal input 

A t = A in sm(2nt/T) (31) 



dv 


1 


dt 


"7 


dw 


€ 


dt 


= 1 



Springer 



Page 12 of 41 



H.G. Rotstein 



has the same frequency (1 cycle per T units of time) for all values of the input fre- 
quency /. The latter affects the speed of the trajectories in the phase plane (i.e., the 
speed of the system's response) without affecting the direction of the underlying vec- 
tor field. Here we use T = 1000. 

The v- and w-nullclines for the unforced system (Aj n = A t = 0) are given by 
N v (v) = — v and N w (v) = av, respectively. The stability properties of this au- 
tonomous system have been discussed in Sect. 3.1.1. For A- m > 0, system (29)-(30) 
is no longer two-dimensional. In our analysis, we will think of the projection of the 
zero-level curve for (29) for each value of t as the f-nullcline N v (v) for the au- 
tonomous system forced by the sinusoidal input A t : 

N Vtt (v) = -v + A t . (32) 

For t = 0 (or any multiple of t = 500), N v j coincides N v = — v. As t increases N Vft 
"moves" cyclically, generating lines parallel to N v = — v in between the lines 

N+(v) = -v + A ia and N~(v) = -v - A- m . (33) 

The moving f-nullcline reaches these two lines at t = 250 and t = 750, respectively. 
As N v ^ t moves so does its intersection with the w-nullcline N w generating a "moving 
fixed point" 

/ A t a A t \ 
(v t9 w t ) = ( T ^ 9 —^) 9 (34) 

which oscillates with frequency 1 between the endpoints (1)250 > ^250) (^250 = ^in) 
and (1)750, wi5o) (A750 = -Am), reaching the origin three times within a cycle at 
t = 0,t = 500, and t = 1000. 

This moving fixed point together with the moving nullclines organize the dynam- 
ics of the forced system. Specifically, the stable fixed point (either a node or a focus) 
acts as a moving target for trajectories, which track its motion by evolving with a 
/-dependent speed. 

The dynamics of the forced system (29)-(30) results from the interaction between 
the vector field of the autonomous system (A[ n = A t = 0) and the forcing term (A t ) 
that causes the cyclic motion of both the i;-nullcline and the fixed point (34) in ad- 
dition to a cyclic change in the direction field. For each point in the phase plane 
the value of the direction field is independent of /. However, the trajectories' f- 
dependent speed causes them to sweep different distances in a given unit of time. 
Since they reach different points, subject to different values (magnitude and direc- 
tion) of the vector field, trajectories describe different curves for different values of /. 
The resonant frequency, if it exists, is the value of the input frequency / for which 
the corresponding limit cycle trajectory has the maximal amplitude in the v -direction, 
provided this amplitude is larger than the instantaneous amplitude corresponding to 
the effective resistance Z(0). 

3.3 Transient Dynamics for Instantaneously Forced Systems 

As a first step in our investigation we describe the response of system (29)-(30) to 
instantaneous constant inputs (with no associated time scales) in order to understand 
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how the system's intrinsic time scales depend on the parameters a and e, and how 
these parameters contribute to the amplification or attenuation of the voltage response 
to these instantaneous changes. We will use this knowledge to explain the dynamics 
of the forced system when the inputs do have associated times scales. 

For the purpose of this part we will consider constant values of A t in system (29)- 
(30). Without loss of generality, we assume that (i) A t instantaneously changes from 
A t = 0 to either A t = 1 or A t = — 1, (ii) this change occurs at t = 0, and (iii) / = 1. 

3.3.1 Transient Dynamics and Effective Time Scales 

An instantaneous change in A t causes an instantaneous displacement of the v- 
nullcline on the phase plane, and consequently an instantaneous change in the loca- 
tion of the fixed point. The vector field changes accordingly, and can also be thought 
of being instantaneously displaced. 

The stability properties of the fixed point, however, remain unchanged since the 
system (29)-(30) is linear and, therefore, its associated Jacobian matrix is indepen- 
dent of A t . 

Consequently, for different values of A t , trajectories at a given initial point have 
the same long term behavior as they approach the corresponding fixed points, but 
they have different initial transient behaviors (Fig. 4). For the purpose of this paper, 
we are particularly interested in the initial transient behavior caused by changes in 
A t . As we explain in more detail below, this transient behavior rather than the long 
term behavior is the one that plays an important role in determining the amplitudes 
of the limit cycle trajectories for the time-dependent forced systems. 

The trajectories' transient behavior relevant to this paper corresponds to the time 
elapsed from the initial point (red dots in Fig. 4) until they cross the corresponding 
f-nullciine (dashed-red line), reaching their maximum value for v. The effect of dif- 
ferent values of A t on the behavior of trajectories initially at the same point can be 
understood by looking at the horizontal and vertical components of the vector field 

D v = — v — w + A t and D w = e(av — w), (35) 

respectively. As A t increases, D v increases and the motion in the horizontal direction 
strengthens relative to the motion in the vertical direction. Consequently, the larger A t 
the "more horizontal" the trajectories' transient direction of motion. From a different 
perspective, trajectories starting at the same initial point in the "standard" coordinate 
system whose origin is (0, 0), start at different initial points in translated coordinate 
systems whose origins are the fixed points (v t , w t ) determined by A t through (34). 
The larger A t , the larger the distance between the trajectory's initial point and the 
corresponding u-nullcline, and so they are less affected by it and they can transiently 
move in more horizontal directions. 

Figure 4 also illustrates the roles of a and € in determining both the effective time 
scales for system (29)-(30) and the amplitude of the voltage response (maximum 
value of v for any given trajectory) to constant perturbations (A t ). For e = 0.01 (top 
panels) the time scales are well separated and the system is fast-slow. Trajectories 
move first fast in almost horizontal directions, then they slow down as they approach 
the u-nullcline, and finally they reverse direction and move towards the fixed point 
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Fig. 4 Phase-plane diagrams for the autonomous linear system (29)-(30) for various representative values 
of a, 6 and At. a a = 1. Top row: e = 0.01. Middle row: e =0.1. Bottom row: e = 1. b a = 6. Top row: 
6 = 0.01. Middle row: e = OA. Bottom row: e = 1. c a = —2. Top row: e = —0.01. Middle row: e = —0.1. 
Bottom row: e = —0.5. Each parce/ shows superimposed phase-planes diagrams for three different con- 
stant values of A t (= 0, 1, — 1) generating three u-nullclines {solid-red for A t = 0, dashed-red for A ? = 1 
and At = — 1) and three fixed points doto on the intersections between the red and greerc /mes). The 
w;-nullcline (green line) is common to all values of A t . Solid-red line: u-nullclines for A t = 0. Dashed-red 
lines: u-nullclines for A t = 1 (above) and A t = —1 (below). Red dots at (0, —1) and (0, —0.5): represen- 
tative initial conditions. Solid-blue lines: trajectories initially located at these initial points. Each trajectory 
emerging from the red dots corresponds to a different value of A t and converges to the corresponding fixed 
point. The fixed points in panels a-top, a-middle and b-top are stable nodes and the fixed points in panels 
a-bottom, b-middle and b-bottom are stable foci 



in close vicinities of the displaced u-nullcline. This effect is more pronounced for 
smaller values of a (compare top panels a and b) because of the product a 6 in (30), 
which affects the effective time scales. As a decreases so does the slope of the v- 
nullcline, thus "pressing" trajectories that become more "compacted". As 6 increases 
(from the top to the bottom panels), the separation of time scales decreases, trajecto- 
ries become more "rounded", and reach lower maximal values of v. 

An important observation is that, qualitatively, the initial transient behavior of 
trajectories in Fig. 4 is not dependent on whether the fixed point is a stable node or 
focus but rather on the magnitude of the differences between the corresponding values 
of a and €. Specifically, trajectories with close enough values of both a and € will 
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Fig. 5 Voltage traces for the autonomous linear system (29)-(30) for various representative values 
of a, 6 and At = 1. a a = 1. The maximum times are ^max, 6=0.01 — 4.31, £ m ax,6=0.1 — 2.38, and 
^max,e=l = 1.12. b 6 = 0.1. The maximum times are t miiX a= \ = 2.38 and t mSLX ,a=6 — 1-44 



display qualitatively similar transient behavior regardless of the possible differences 
in the stability of the corresponding fixed points (Fig. 5). 

3.3.2 Amplification of the Voltage Response to Instantaneous Constant Inputs 

The voltage response is amplified by a given parameter if v is able to reach a higher 
value when this parameter changes. Figures 4 and 5 demonstrate that the amplifica- 
tion of the voltage response increases as the time- scale separation increases (e de- 
creases) (compare top, middle and bottom panels for increasing values of € in Fig. 4 
and see the traces in Fig. 5) and as a decreases (compare Figs. 4a and 4b and see the 
traces in Fig. 5). These results reflect the roles of the resonant and amplifying currents 
in determining the amplification of the voltage response through the dimensionless 
parameters a and e given by (12). For example, an increase in time constant T\ causes 
a decrease in € (time-scale separation), leading to an amplification of the voltage re- 
sponse. The explanation of the effects of other biophysical parameters (e.g., #l and 
gi) through the dimensionless parameter a is less straightforward since a has been 
used as a voltage scale. A thorough study of these effects has been carried in [29]. 

Figure 4 also demonstrates that voltage responses are more amplified for negative 
values of a and € than for positive ones (compare panels a and c). This is consis- 
tent with the fact that negative values of a and € are obtained for negative values 
of the effective conductance which in turn reflects the presence of a strong, fast 
amplifying current (h) in the biophysical model. 

3.4 Voltage Amplitude Response to Sinusoidal Inputs at Different Frequencies: 
Resonance and Phase-Resonance 

Here we build up on the ideas discussed in Sect. 3.3 to investigate the mechanism 
of selection of the resonant and phase-resonant frequencies in the two-dimensional 
system (29)-(30) in response to sinusoidal inputs of the form A t (31). Without loss 
of generality we consider the canonical case A- m = 1. With this choice, for each 
value of / the impedance function (18) coincides with the maximum value of v\ 
i.e., = Z(f)> 

Ast increases, At changes and the f-nullcline N v t (32) moves cyclically, parallel 
to itself between the lines N*(v) and N~(v) (33) (dashed-red lines in Fig. 4) with 
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frequency equal to 1 . The fixed point (34) moves cyclically with the same frequency 
between the two corresponding extreme values. The input amplitude is coded by the 
distance between the moving fixed point and the origin. 

Trajectories respond to changes in A t by evolving according to the dynamics dic- 
tated by the underlying vector field (35). The trajectories' speed v is given by the 
product of the magnitude of the vector field and f~ l 



At any given point, the underlying vector field is independent of /. However, the 
/-dependent speed (36) causes the resulting limit cycle trajectory to be /-dependent 
as explained at the end of Sect. 3.2. This dependence is reflected primarily on the 
shape and orientation of the limit cycle as we illustrate this in Fig. 6a for a system 
that exhibits resonance (see Fig. 6b) and phase-resonance (see Fig. 6c). For small 
values of / limit cycles are elliptic and very narrow, with their major axis lying on 
a vicinity w-nullcline. As / increases, the limit cycle first widens and rotates, then 
it narrows as it continues to rotate until the major axis becomes horizontal (on the 
i;-axis), and finally the limit cycle shrinks until it collapses to a point as / — ► oo. The 
system exhibits resonance because the maximal value of v for / ~ 65 is larger than 
the maximal value of v for / —> 0. 

3.4. 1 Envelope -Plane Diagrams: Resonant and Phase -Re sonant Responses 

Figure 6d shows a graph containing the curve generated by the points with maximum 
values of v on the limit cycles for continuously changing values of / e [0, oo) and 
the parameter values in Fig. 6a. We refer to this curve as the upper envelope curve 
of the voltage response. (The lower envelope curve, not shown, is determined by 
the minimum values of v on the /-dependent limit cycles, and is symmetric to the 
upper envelope curve with respect to the f-nullcline and w-nullclines.) We refer to 
the diagrams containing the envelope curves together with the v- and w-nullclines 
(green, solid-red, and dashed-red), as in Fig. 6d, as the envelope -plane diagrams. 

Envelope-plane diagrams contain geometric and dynamic information about a sys- 
tem's frequency response to oscillatory inputs, and they are the frequency analogs to 
phase-plane diagrams. Trajectories in the envelope-plane diagrams (upper and lower 
envelopes) are curves parameterized by the input frequency as trajectories in the 
phase planes are curves parametrized by time. Neither / nor t are explicit in the 
corresponding diagrams. The red-dashed lines quantify the maximal input amplitude. 
For / — >► 0, the envelope curve is at the intersection between the upper dashed-red 
line and the w-nullcline. For f —> oo, the envelope curve is at the origin (limit 
cycle shrinking to a point). 

The cusp ("horizontal peak") in the envelope curve in Fig. 6d corresponds to the 
peak in the impedance profile (Fig. 6b), and hence to the resonant frequency. At this 
cusp the voltage response is the largest across input frequencies, and larger than the 
response for /—►(): v max = (1 + a) -1 . 

The point in the envelope curve tangent to the upper dashed-red line corresponds 
to the phase-resonant frequency for reasons that we will clarify later in the paper. 



v = 



v — w + A t ) 2 + € 2 (av — w) 2 



(36) 
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Fig. 6 Dynamics of the sinusoidally forced linear system (29)-(30) for a = 1, e = 0.01. a Projections 
of the phase-plane diagrams on the v-w plane for various representative values of the input frequency /. 
The fixed solid-red line and the solid-green line represent the v- and u;-nullclines for the unforced sys- 
tem (A t = 0), respectively. The dashed-red lines represent the u-nullclines for A t = 1 (above-right) and 
At = — 1 (below-left). The solid-blue lines represent the trajectories of the forced system for a single pe- 
riod (T = 1000). b Impedance profile, c Phase profile, d Envelope-plane diagram. The solid-blue line 
represents the envelope curve: Each point on this curve is the maximum point on the limit cycle response 
to sinusoidal inputs parametrized by the input frequency / which increases from f = 0 (blue-square at 
the intersection between upper dashed-red and green curves) to / -> oo (blue dot at the origin). (The 
v -coordinates of the envelope curve are the impedance function Z, since A[ n = 1.) Solid-red and -green 
lines: v- and u;-nullclines for the unforced system (A t = 0). Dashed-red lines: u-nullclines for A t = 1 
(above) and A t = —1 (below) 

Below, we first discuss the voltage response mechanisms for the two limiting cases 
(/ -> 0 and f oo) and then elaborate on the dynamics for intermediate values of 
/, in particular / res and / p h as . 

3.4.2 Low Input Frequencies Generate Quasi-one-dimensional Dynamics Along the 
w-Nullcline 

For values of / 1 both equations in system (29)-(30) are very fast (dv/dt — ► oo 
and dw/dt -> oo), and then the limit cycle trajectory tracks the motion of the fixed 
point almost instantaneously (in a quasi- steady- state fashion). In the limit /—►(), the 
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Fig. 7 Snapshots of the moving phase-plane diagram for the linear system (29)-(30) for a = 1 , e = 0. 1 , 
f = 1 and representative values of t within one cycle (T = 1000). The solid-red lines represent the mov- 
ing u-nullcline. The dashed-red lines represent the u-nullclines for A t = 0 (middle), A t = 1 (top), and 
A t = — 1 (bottom). The solid-green line represents the (fixed) w-nullcline. The fixed solid-red line and 
the solid-green line represent the v- and u;-nullclines for the unforced system (A t = 0), respectively. The 
dashed-red lines represent the u-nullclines for A t = 1 (above), A t = 0 (middle), and A t = — 1 (below). As 
t increases, the red nullcline moves cyclically between the two dashed-red lines. The blue dot indicates the 
initial location of the trajectory (t = 0 = 1000) 



limit cycle trajectory moves cyclically along the w -nullcline in between the points 
(34) with A t = ±Ai n = ±1. From (34), fmaxMin = (1 + a) -1 , which coincides with 
Z(0) in (21). In Fig. 7 we present snapshots of the evolution of the limit cycle trajec- 
tory for / = 1 (corresponding to Fig. 6a, top-left panel). 

3.4.3 High Input Frequencies Generate Quasi- one -dimensional Dynamics Along 
the Horizontal v-Axis 

For large enough values of / ^> 1, both (29) and (30) are very slow (dv/dt -> 0 
and dw/dt -> 0), and then the limit cycle trajectory evolves with a very low speed 
according to (36). Relative to the trajectory's speed, the v -nullcline N v j moves very 
fast. The time it takes the trajectory to cover a very small distance while tracking 
the motion of the v -nullcline is enough for the latter to reach its maximum level 
reverse direction, and intersect the trajectory on its "way back". This causes 
the trajectory to reverse direction. As a result, the limit cycle trajectory is not able 
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to cover any significant distance away from the origin, and so the amplitude of the 
limit cycle is small as compared to other values of / (Fig. 6a, bottom-right panel). In 
the limit f —> oo, the amplitude of the limit cycle trajectory is zero (the limit cycle 
shrinks to the origin), which coincides with lim/-^ Z(f) = 0. 

3.4.4 Voltage Response Amplification at the Resonant Frequency Band 

For small and large enough values of / the corresponding limit cycle trajectories 
are constrained to move along quasi-one-dimensional directions (w-nullcline and v- 
axis, respectively). For intermediate values of / there is a transition in the shapes 
of the limit cycle trajectories between these two limit cases. Specifically, limit cycle 
trajectories are neither too fast nor too slow, and then, while they are "left behind" 
by the moving u-nullcline, they can take advantage of the two-dimensional vector 
field without being constrained to move in quasi-one-dimensional directions. For the 
appropriate values of a and € this degree of freedom allows limit cycle trajectories to 
reach values of v m2LX (f) larger than f m ax(0), and so to exhibit resonance. 

In Fig. 8 we show a sequence of snapshots of the evolution of the limit cycle trajec- 
tory for / = / res = 65 (Fig. 6a, bottom-left). The initial snapshot (t = 0) corresponds 
to a point (blue dot) on the limit cycle trajectory and A t = 0. As t increases the v- 
nullcline moves and the limit cycle trajectory tracks its motion. The limit cycle' shape 
and amplitude result from the combined effect of A t , f, and the model parameters 
(a and e). 

In Fig. 8 we are using a relatively small value of € (e = 0. 1). While, small values of 
€ are enough to account for the almost horizontal transient trajectories in autonomous 
systems (Fig. 4), they are not enough to account for the almost horizontal directions 
of motion for limit cycles trajectories in forced systems such as in Fig. 8. In fact, for 
the same value of € but a different value of / (/ = 1 instead of / = 65) the limit 
cycle trajectory moves in the (oblique) direction of the w-nullcline (Fig. 7). 

To better understand how A t , f, € and a affect the direction of motion of the 
limit cycle trajectory in Fig. 8, it is useful to go back and look at the initial transient 
segment of the evolution of trajectories for the autonomous systems presented in 
Fig. 4 (from the initial, red point until they cross the f-nullcline). Figure 4a (middle) 
corresponds to the same parameter values (a = 1 and € =0.1) as in Fig. 8. The 
increase in the constant input from A t = 0 to A t = 1 causes both an increase in the 
distance between the tip of the trajectory and the u-nullcline (from the solid line to 
the dashed-red line) and a strengthening of the horizontal component of the vector 
field D v (35) relative to its vertical component D w . As a result, during the initial 
transient interval, the larger A t > 0, the more horizontal is the trajectory's direction 
of motion. 

A similar, but dynamic effect is partially responsible for the determination of the 
direction of motion in Fig. 8. Trajectories move faster in more horizontal directions 
as A t increases during the ascending phase (t = 0 to t = 250) because the contin- 
uous increase in A t causes a continuous increase in D v while D w remains almost 
unchanged. However, differently from the static case, in the dynamic case there is 
an opposite, dynamic effect caused by the decreasing distance between the tip of the 
limit cycle trajectory and the moving u-nullcline N V) t since D v decreases as the limit 
cycle trajectory approaches the f-nullcline. 
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Fig. 8 Snapshots of the moving phase-plane diagram for the linear system (29)-(30) for a = 1 , e = 0. 1 , 
/ = 65 and representative values of t within one cycle (T = 1000). The solid-red lines represent the 
moving u -nullcline. The dashed-red lines represent the u-nullclines for A t = 0 (middle), A t = 1 (top), 
and At = — 1 (bottom). The solid-green line represents the (fixed) w-nullcline. The fixed solid-red line and 
the solid-green line represent the v- and u;-nullclines for the unforced system (A t = 0), respectively. The 
dashed-red lines represent the u-nullclines for A t = 1 (above), A t = 0 (middle), and A t = — 1 (below). As 
t increases, the red nullcline moves cyclically between the two dashed-red lines. The blue dot indicates the 
initial location of the trajectory (t = 0 = 1000) 



In the autonomous case, the nullclines are fixed, so the distance between the tra- 
jectory and N Vit depends only on the trajectory's speed (36). In the forced case, on 
the other hand, the trajectory and the v -nullcline N v j are simultaneously moving, 
with different speeds. The resonance frequency / res is the input frequency for which 
their relative speed is optimal in the sense that it allows the limit cycle trajectory to 
move in a direction that maximizes (over the range of input frequencies /) the dis- 
tance (in the v direction) the limit cycle trajectory can cover before intersecting the 
moving v -nullcline, at which time they are forced to reverse direction. This causes 
a maximization (over the range of input frequencies /) of the maximum value t> max 
on the limit cycle trajectories. In Fig. 8 (/ ~ /res), this intersection occurs for t> max 
(~ 1). The remainder of the limit cycle can be explained using similar ideas. 

For values of / < / res , the limit cycle trajectory moves faster than for / = / res , 
and so the direction of motion is less horizontal causing this cycle trajectory to inter- 
sect N V) t at a higher point; i.e., for a lower value of f max (e.g., Fig. 6a, / = 48). For 
values of / > / res , the limit cycle trajectory moves slower than for / = / res . The rela- 
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tive speed between this trajectory and N VJ is large enough to allow the limit trajectory 
to move in an almost horizontal direction, but, due to the limit cycle trajectory's lower 
speed, N Vi t "returns" before this trajectory covers a large enough distance. Thus, the 
intersection between the trajectory and N v , t occurs for lower values of t> max (e.g., 
Fig. 6a, / = 200). 

3.4.5 Resonance and Intrinsic Oscillations Are Generated by Related, but not 
Identical Mechanisms 

Intrinsic STOs and subthreshold resonance have been proposed to result from the 
same underlying mechanism [12]. However, in linear systems, resonance may occur 
in the absence of intrinsic oscillations and vice versa [11] (see Fig. 3a). Even for 
parameter values for which linear systems exhibit both resonance and intrinsic oscil- 
lations, the resonant and natural frequencies do not necessarily coincide (Fig. 3b). 

The phase-plane analysis discussed above demonstrates that the resonant proper- 
ties are not expected to be directly linked to the stability properties of the underlying 
autonomous system. The time scales corresponding to intrinsic oscillations, if they 
exist, are determined by the stability properties of the stable foci, which reflect the 
long term behavior of trajectories. In contrast, the time scales corresponding to the 
resonant frequency do not involve the stability properties of the fixed points of the un- 
derlying autonomous system (nodes or foci) but rather the initial transient behavior 
of trajectories as described above. This initial transient behavior is the link between 
the autonomous and forced systems. This transient dynamics is associated with ei- 
ther the onset of intrinsic oscillations (foci) or a "sag" (node) typically observed in 
the response of h-currents to constant inputs [11]. Importantly, as we demonstrated 
in Sect. 3.3, autonomous systems may display qualitatively similar transient behavior 
for nodes and foci, even though the long time behavior of the corresponding trajecto- 
ries is qualitatively different (see Fig. 4b). 

3.4.6 Phase Advance, Phase Delay and Phase -Resonance 

In Fig. 3b we showed that the resonant (/res) an d phase-resonant (/ p has) frequencies 
do not necessarily coincide, and that resonance may occur in the absence of phase- 
resonance. For the parameters values in Fig. 6, / res = 65 and / p has = 48. Typical 
phase profiles show "phase advance" for / < / p h as and "phase delay" for / > / p h as 
(see Fig. 6c). 

Geometrically, for the output and input to be synchronized in phase, the intersec- 
tion between the limit cycle trajectory and the moving u-nullcline must occur when 
the i>-nullcline reaches its highest level A t = A[ n (when the solid-red line in the mov- 
ing phase-plane diagrams reaches the dashed-red line at t = 250). The phase is ad- 
vanced when the limit cycle trajectory "peaks" before the input does, which requires 
the trajectory to evolve fast enough as compared to the speed of the f-nullcline. Con- 
versely, the phase is delayed when the limit cycle trajectory intersects the u-nullcline 
"on its way back" (descending phase), after the u-nullcline has reached its high- 
est level. For / = / res = 65 (Fig. 8), the limit cycle trajectory is still behind the 
f-nullcline at the time this u-nullcline reaches its highest level (t = 250), and the 
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Fig. 9 Snapshots of the moving phase-plane diagram for the linear system (29)-(30) for a = 1 , e = 0. 1 , 
/ = 48 and representative values of t within one cycle (T = 1000). The solid-red lines represent the 
moving u -nullcline. The dashed-red lines represent the u-nullclines for A t = 0 (middle), A t = 1 (top), 
and At = — 1 (bottom). The solid-green line represents the (fixed) w-nullcline. The fixed solid-red line and 
the solid-green line represent the v- and u;-nullclines for the unforced system (A t = 0), respectively. The 
dashed-red lines represent the u-nullclines for A t = 1 (above), A t = 0 (middle), and A t = — 1 (below). As 
t increases, the red nullcline moves cyclically between the two dashed-red lines. The blue dot indicates the 
initial location of the trajectory (t = 0 = 1000) 



phase is therefore delayed (the intersection occurs at a time in between the top-right 
and middle-left panels, t = 250 and t = 375, respectively). 

For / = /phas = 48 (Fig. 9), on the other hand, the limit cycle trajectory intersects 
the v -nullcline when the latter is at its highest level (t = 250). Since the point of 
intersection corresponds to v = L> max (on the limit cycle trajectory), the correspond- 
ing point on the envelope curve "touches" the top dash-red line. For all other fre- 
quencies, the intersection between limit cycle trajectories and the moving v -nullcline 
occurs away from the top dashed-red line, therefore the phase-resonant frequency 
corresponds to the point in the envelope curve that is tangent to the top dashed-red 
line. This tangent point is not necessarily the cusp that corresponds to the resonance 
frequency / res . 

Specifically, the phase-resonance phenomenon depends on the ability of the limit 
cycle trajectory for / = / p h as to track the v -nullcline fast enough so to reach t> max 
before the v -nullcline reaches its highest level. But this does not preclude limit cycle 
trajectories for frequencies / ^ / p has to reach higher values of v due to the different 
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directions of motion generated by these values of /. In fact, values of / > / p has 
cause limit cycle trajectories to evolve slower than for / = / p has, and so (within 
some range of values of /) they can move along more horizontal directions, and thus 
reach higher values of v before intersecting the "returning" u-nullcline. This is the 
case for / = / res in Fig. 8. 

3.5 Effects of a and € on Resonance, Phase-Resonance, and Other Attributes of the 
Impedance and Phase Profiles 

Geometrically, the value of a determines one of the boundaries of the region (triangle) 
in the envelope-plane diagram where the envelope curves live (Fig. 6d). The other two 
boundaries are the u-nullcline for A t = 1 (dashed-red line) and the u-axis. The value 
of a also determines Z(0) = f m ax(0)Min = (1 + a) -1 , which is the v -coordinate 
of the intersection between the w-nullcline and the u-nullcline for A t = 1. As we 
mentioned earlier, for resonance to occur, v mdiX (f) > v maK (0) (on the envelope curve) 
for some range of values of / around / res . 

The value of € plays a key role in determining the direction of motion of limit cycle 
trajectories analogous to its effect on the initial transient behavior of trajectories for 
the autonomous systems discussed in the context of Fig. 4. For fixed values of a, the 
smaller € the more horizontal (and less rounded) is the direction of motion of the 
initial portion of the trajectories, and the larger the value of v these trajectories reach; 
i.e., the voltage response is amplified. 

From the biophysical point of view, a and € convey information about the effec- 
tive conductances g\ and the capacitance C, and the time constant x\ (gating 
variable x\) through (12). In [29] we use numerical simulations to conduct a through 
analysis of the effects of resonant and amplifying biophysical conductances on the 
determination of / res , / pri as and other attributes of the impedance and phase profiles. 

Here we use the envelope-plain diagrams developed in Sect. 3.4 to explain how 
the parameters a and € affect the resonant and phase-resonant properties of (14)- 
(15), and how they contribute to the amplification of the voltage response. In Figs. 10 
and 11 we examine the effects of changes in € and a, respectively. In Figs. 12 and 13 
we investigate the additional amplifying effects for negative values of both a and €. 

3.5.1 An Increase in the Time-Scale Separation (Decreasing Values ofe) Causes an 
Amplification of the Voltage Amplitude Response 

In Fig. 10 we compare the voltage responses for various representative values of € 
and a = 1. As € decreases, both Z max and <2z increase, / res and / pri as decrease, and 
the neuron becomes more selective (A\/2 decreases) (Fig. 10b). These differences are 
captured by the shapes of the envelope curves, which are sharper the smaller the €. 
For € = 0.01 (panel a, left) the system is fast- slow causing trajectories to move fast 
along horizontal fibers as compared to larger values of € (middle and right panels), 
and thus the limit cycle trajectory corresponding to / res reaches a higher value of 
fmax (the voltage response is more amplified). The envelope-plane diagrams predict 
that / r es and / p has are very close since the envelope curve is tangent to the dashed-red 
line almost at the peak (panel a). 
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As € increases, the tangent point and the peak separate (panels b and c), and hence 
the difference between / res and / p has increases. For € = 1 (panel c), the tangent point 
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Fig. 12 Dynamics of the sinusoidally forced linear system (29)-(30) for a = —2, e = —0.5. a Projections 
of the phase-plane diagrams on the v-w plane for various representative values of the input frequency /. 
The fixed solid-red line and the solid-green line represent the v- and u;-nullclines for the unforced sys- 
tem (A t = 0), respectively. The dashed-red lines represent the u-nullclines for A t = 1 (above-right) and 
At = — 1 (below-left). The solid-blue lines represent the trajectories of the forced system for a single pe- 
riod (T = 1000). b Impedance profile, c Phase profile, d Envelope-plane diagram. The solid-blue line 
represents the envelope curve: Each point on this curve is the maximum point on the limit cycle response 
to sinusoidal inputs parametrized by the input frequency / which increases from f = 0 (blue-square at 
the intersection between upper dashed-red and green curves) to / -> oo (blue dot at the origin). (The 
v -coordinates of the envelope curve are the impedance function Z, since A[ n = 1.) Solid-red and -green 
lines: v- and u;-nullclines for the unforced system (At = 0). Dashed-red lines: u-nullclines for A t = 1 
(above) and A t = — 1 (below) 

coincides with the intersection between the green and dashed-red lines, and hence the 
system exhibits no phase-resonance (/ p has = 0), while it still exhibits resonance. 

In the limit of € —> 0, the system becomes quasi-one-dimensional and the dynam- 
ics occurs almost exclusively along the u-axis for almost all values of /, except a 
small range close to / = 0. For € = 0 the dynamics is fully one-dimensional, thus 
generating a low-pass filter response (no resonance). For / — ► 0 (and € = 0) in (29) 
and (30) the maximum value of v approaches the intersection between the dotted- 
red line and the u-axis in the envelope-plane diagram. The maximum value of v 
(along the u-axis) decreases as / increases. As soon as € > 0, the dynamics becomes 
two-dimensional. The envelope curves unfold in the two-dimensional space and ac- 
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Fig. 13 Effects of changes in e on the resonant properties of the linear system (29)-(30) for a = —2. 
a Envelope-plane diagrams for e = —0.01 (left), e = —0.1 (middle), and e = —0.5 (right), b Impedance 
(Z) profiles, c Phase (</>) profiles 



quire "triangular-like" shapes similar to these shown in Fig. 10a (left), but peakier the 
smaller €. Regardless of how small is €, there is always a range of values of / < e 
such that as / —> 0, the limit cycle trajectory evolves fast, "almost along" the w- 
nullcline (green line) following the motion of the fixed point. Therefore, the point on 
the envelope curve for / = 0 lies on the intersection of the dotted-red line and the w- 
nullcline. The remaining of the envelope curve will be as in Fig. 10a, but the distance 
between the envelope curve and the dotted-red line will be larger the smaller € . This 
distance is determined by the minimum phase, which increases in absolute value as € 
decreases. 

3.5.2 Increasing Values of a Generate Resonance and Phase-Resonance and Cause 
an Amplification of the Voltage Response 

In Fig. 1 1 we compare the voltage responses for various representative values of a 
and 6 = 1. The value of a determines the slope of the w-nullcline, which is steeper 
for larger values of a, and f m ax(0) = (1 + a) -1 . The smaller slope of the w-nullcline 
for lower values of a (right panel) reduces the trajectories' freedom of motion by 
constraining them to evolve in close vicinities of both nullclines (the w-nullcline and 
the moving u-nullcline). The consequent decrease in speed prevents the limit cycle 
trajectory from reaching large enough values of v mdx - This together with the fact that 
fmax(0) increases with a prevents the system from exhibiting resonance. 

For larger values of a (left and middle panels), v max (0) is smaller. Although an 
increase in a increases the horizontal component of the vector field (through the 
product ae) and hence reduces the effective time-scale separation, this is not enough 
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to prevent limit cycle trajectories from moving beyond v mSLX (Q), and the system is 
able to exhibit resonance. 

The envelope-plane diagrams predict the generation of phase-resonance as a in- 
creases above some critical value (in between these for a = 2 and a =0.5). For 
a = 0.2 and a = 0.5 (middle and right panels) the envelope curve is tangent to the 
dashed-red at its intersection with the w-nullcline, and hence the system exhibits no 
phase-resonance. In contrast, for a = 2 (left panel), the tangent point occurs in be- 
tween this intersection and the peak of the envelope curve, and thus the system does 
exhibit phase-resonance. 

3.5.3 Negative Values of Both a and € Amplify the Voltage Response 

From (12), negative values of both a and € reflect the presence of a strong enough 
amplifying current that makes the effective conductance #l < 0 and a resonant cur- 
rent (g\ > 0). For the underlying autonomous system to be stable the values of € are 
constrained to be in the region € > — 1 (see Fig. 3a). We present a representative ex- 
ample in Fig. 12 (or = —2 and € = — 0.5). Comparison with the examples presented 
in Figs. 10 and 11 demonstrates that the voltage response is more amplified for neg- 
ative than for positive values of both a and € . This is true even for negative values 
of € that are not very small in absolute value (time-scale separation no necessarily 
small). 

Geometrically, negative values of a cause the slow of the w-nullcline to be nega- 
tive (Fig. 12a) and affect the direction of motion of trajectories. The dynamics of the 
forced system for very low and very large values of the input frequency / is similar 
to the cases discussed above for positive values of both a and € (Fig. 12 for / = 3 and 
/ = 400, respectively). For intermediate values of /, the amplification of the voltage 
response reflects the interaction of the oscillatory inputs with an intrinsic vector field 
(for the autonomous system) that is correlated with the amplification of the voltage 
responses to instantaneous inputs discussed in the context of Fig. 4c for negative as 
compared to positive values of both a and € (Figs. 4a and 4b). 

There are three main differences between the two cases (positive and negative 
values of both a and €). First, the phase-resonant frequency / p h as = 138 is larger than 
the resonant frequency / res = 108 for negative parameter values (Fig. 12b), while 
/phas < /res for positive parameter values. This difference is also captured by the 
envelope-plane diagram in Fig. 12d. Second, the voltage response is amplified by 
decreasing, rather than increasing levels of the time-scale separation (given by |e|), 
for negative parameter values (Figs. 13a and 13b). Finally, the phase profiles are 
initially at 0 = — re and increase towards 0 = tt/2 (Fig. 13c) for negative values of 
both a and € (compare with Figs. 10 and 11), while 0 never decreases below —tt/2 
for positive values of these two parameters. 

3.6 Extension of the Envelope-Plane Diagram Approach to Simple Nonlinear 
Systems 

Linearization around resting potential produces a good approximation to the 
impedance and phase profiles for weakly forced linear systems. However, signifi- 
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cant departures from the linear prediction are expected for larger input amplitudes in 
the presence of strong nonlinearities, in particular when they are close to the voltage 
threshold for spike generation. The question arises of how nonlinearities affect the 
voltage response and, in particular, how the resonant and phase-resonant frequencies 
depend on the type of these nonlinearities and the time- scale separation between the 
participating variables. 

Linear responses to sinusoidal inputs are characterized by (i) the coincidence of 
the input and output frequencies, (ii) the proportionality between the output and the 
input signals that renders the impedance function independent of the input ampli- 
tude Aj n , and (iii) the symmetry between positive and negative values of the voltage 
response. In terms of the envelope-plane diagrams, the information about the volt- 
age response is captured by the upper envelope curve for A- m = 1 . Here we illustrate 
how the ideas developed in the previous sections can be extended to the investigation 
of the mechanisms underlying the generation of resonance and phase-resonance in 
nonlinear systems. 

We use simple piecewise-linear (PWL) systems of the form 

h v (v) - w + A in sm(fit), (37) 

€[h w (v)-w], (38) 

where the functions h v and h w are continuous PWL functions with two linear pieces. 
We consider two representative cases! (i) h v is PWL and h w is linear, and (ii) h v 
is linear and h w is PWL. In both cases, the nullclines intersect at the origin. By 
construction, in a vicinity of the fixed point (v, w) = (0, 0) both h v and h w are linear, 
and system (37)-(38) has the form (14)-(15). 

3.6.1 The Voltage Response Is Amplified by Nonlinearities in the Voltage Equation 
Here we consider system (37)-(38) with 

Irjv ift><0.8, 
and h w (v)=av (39) 
r\ r v \tv> 0.8 

with T) = — 1, r\ r = —0.4 and a = 1. Our results are presented in Fig. 14. Figure 14a 
shows the envelope-plane diagrams (including both the upper and lower envelopes) 
for representative values of a and e. The f-nullcline for the unperturbed system 
(solid-red lines) breaks at v = 0.8. For the forced system, the location of the v- 
nullcline changes with t and so does the breaking point. The impedance and phase 
profiles for the parameters in Figs. 15a (€ = 0.01) and 15b (e = 0.1) are presented in 
Figs. 15dl and 15d2, respectively. 

For small values of A- m (Figs. 15al and 15bl) the dynamics of the PWL system 
is governed by the underlying linear system and the voltage response is linear. The 
envelope-plane diagrams are symmetric and similar to the ones discussed in previ- 
ous sections. For larger values of A- m the voltage response exhibits an amplification 
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Fig. 14 Voltage response of the PWL system (37)-(38) with h v and h w given by (39) to sinusoidal 
inputs, a-c Envelope-plane diagrams for representative values of A m and e. al A in = 0.8, e = 0.01. 
a2 A in = 1,6 = 0.01. a3 A in = 1.2, e = 0.01. bl A in = 0.8, 6 = 0.1. b2 A in = 1.2, e = 0.1. c A in = 1.2, 
6 = 1. d Impedance profiles for representative values of A m and e. dl e = 0.01. d2 € = 0.1. We used the 
following parameters: a = 1, rj = — 1, rj r = —0.4, v c = 0.8 



that is not captured by the linearization (Figs. 15dl and 15d2). As A- m increases, the 
values of both Z max and Qz increase, the resonant frequency / res decreases, and the 
selectivity increases. Notably, this nonlinear amplification of the voltage response 
is more pronounced the smaller the value of € (larger levels of time- scale separa- 
tion). 

The nonlinear amplification of the voltage response can be understood by com- 
paring the envelope-plane diagrams in Fig. 14a. Geometrically, the nonlinearities are 
reflected by a sudden change in the slope of the f-nullcline. Similarly to the linear 
case, as t progresses, the nonlinear f-nullcline moves cyclically between the two 
dashed-red curves corresponding to ±Ai n . The envelope curves live in the region 
bounded by the w-nullcline, the moving w;-nullcline (dashed-red line) and the u-axis 
(horizontal w = 0 line). As A- m increases, the area of this region increases as the tri- 
angle in panel al becomes a (non-triangular) polygon in Figs. 15a2 and 15a3. These 
changes are accompanied by an increase in the voltage amplitude of the regions that 
allow trajectories to reach higher values of v. 



Springer 



Page 30 of 41 



H.G. Rotstein 




Fig. 15 Voltage response of the PWL system (37)-(38) with h v and h w given by (40) to sinusoidal 
inputs, a A[ n = 0.8. b A[ n = 1. c A[ n = 1.5. We used the following parameters: a = 1, r] = — 1, a r = 0.4, 
v c = 0.5, 6 = 0.01 



The nonlinearities do not affect the dynamics of the limit cycle trajectory for low 
and large values of /, which behaves as explained above for linear systems, but they 
do affect the dynamics of the limit cycle trajectories for intermediate values of /. 
For Ai n = 0.8 (Fig. 15al), the voltage response is quasi-linear since the nonlinearity 
is not present in the triangular region where the envelope curve is located. The limit 
cycle trajectories cannot move beyond the piece of the f-nullcline with slope rj = — 1 . 
For larger values of A- m the nonlinearity moves above the horizontal axis (Figs. 15a2 
and 15a3) during the ascending phase of the oscillatory input, thus changing the shape 
of the envelope-curve region (to a non-triangular polygon). 

3.6.2 The Voltage Response Is Amplified by the Time-Scale Separation in the 
Presence of Nonlinearities in the Voltage Equation 

The nonlinear amplification of the voltage response results from the ability of the 
limit cycle trajectories to move beyond the dashed-red line with slope r] = — 1. 
Clearly, this phenomenon is more pronounced for larger values of A- m . As expected 
from our previous discussion for linear systems, this phenomenon is also more pro- 
nounced for smaller values of € (compare Figs. 14a and 14b) since limit cycle trajec- 
tories move in more horizontal directions for smaller values of € , and they are able to 
reach larger voltage values. As € increases, limit cycle trajectories move in less hori- 
zontal directions and, consequently, they do not take advantage of the extra available 
region of motion in the envelope-plane created by the nonlinear f-nullcline. 

We emphasize that the nonlinearities are present, and they are the same, for all val- 
ues of €, but the underlying vector field causes the limit cycle trajectories to "ignore" 
these nonlinearities and move inside the linear region (Fig. 14d). 

The asymmetry in the system's nonlinearity is reflected in the asymmetry in the 
voltage responses, which is captured by the differences between the upper and lower 
envelope curves in the envelope-plane diagrams. Clearly, the difference between these 
two curves increases as A- m increases. 

We note that information about this asymmetric voltage response is not captured 
by the impedance function. 



4^ Springer 



Journal of Mathematical Neuroscience (2014) 4:11 



Page 31 of 41 



3.6.3 The Voltage Response Is Almost Insensitive to N online arities in the Gating 
Variable Equation 



Here we consider system (37)-(38) with 



h v = rjv and h w (v) = 



av ifu<0.5, 
a r v if v > 0.5 



(40) 



with x] = — 1, a = 1 and a r = 0.4. Envelope-plane diagrams for representative values 
of € are shown in Fig. 15. The nonlinearities are reflected in the sudden change in 
the slope of the w-nullcline. This change slightly affects the dynamics of the limit 
cycle trajectories for low values of /, since they evolve in close vicinities of the w- 
nullcline but they do not significantly affect the dynamics of the limit cycle trajectory 
for larger values of /. Even for large values of A m (not shown), the nonlinearities 
in the w-nullcline do not cause any significant change in the envelope-curve region 
in the envelope-plane, and thus the impedance function does not significantly change 
with increasing values of Am. The resonant frequency does decrease with increasing 
values of A m but significantly less than for the case considered above (not shown). 



4 Discussion 

Subthreshold (membrane potential) resonance has been investigated in a number of 
neuron types both experimentally and theoretically [10, 12-22, 29, 40], and has been 
implicated in the generation of preferred neuronal firing rates [11] and network oscil- 
lations [23, 25, 27, 41-43]. Phase-resonance has also been investigated in neuronal 
systems [11, 14, 29], although to a lesser extent than subthreshold resonance, and it 
has been proposed to play an important role in neuronal synchronization [30]. 

The voltage response to oscillatory inputs in neurons and other electric circuits is 
typically characterized in terms of the impedance and phase profiles [10]. For lin- 
ear systems receiving sinusoidal inputs, these curves can be computed analytically. 
For nonlinear systems or linear systems receiving non-sinusoidal oscillatory inputs, 
analytical calculations are not possible, and these curves are computed numerically. 

Even when analytical expressions for the impedance and phase are available, the 
information they provide about the dynamic mechanisms leading to resonance is lim- 
ited. For instance, the impedance and phase profiles fail to provide the necessary 
insight into the roles played by the interaction between the neuron's intrinsic time 
scales and the time scales associated with the input currents in the generation of reso- 
nance. Furthermore, linearized models fail to capture important nonlinear properties 
of the voltage response such as its nonlinear amplification for significant levels of the 
time-scale separation (small values of e). 

In this paper we have developed and used dynamical system tools to investigate 
the dynamic mechanisms of generation of subthreshold and phase resonance in two- 
dimensional linear and linearized (conductance-based) models, and we have extended 
these tools to include the effect of simple, but not necessarily weak, types of nonlin- 
earities. The mechanistic analysis and the envelope-plane diagrams we developed in 
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this paper provide a framework for the investigation of the preferred frequency re- 
sponses in three-dimensional and nonlinear models. They can also be used as tools 
complementary to both the numerically computed and the experimentally measured 
impedance and phase profiles. More research is needed to adapt these tools to small 
neuronal networks. 

Mechanistic studies of subthreshold resonance have mainly focused on the role 
that resonant and amplifying currents (and their associated gating variables) play in 
shaping the neuronal voltage response [10, 1 1, 29]. As we have shown in [29], the de- 
termination of the attributes and phase profiles as the result of the interaction between 
these two current types is more complex and not as straightforward as previously 
thought [10] (see our discussion of resonant and amplifying currents in Sect. 2.5). 
This is further emphasized by the fact that resonance may occur in the absence of 
intrinsic oscillations and vice versa [11, 29] (see also Fig. 3). 

In this study, we set out to understand the dynamic mechanisms underlying the 
generation of resonance and phase-resonance for a generic class of linearized bio- 
physical models. We were particularly interested in issues that cannot be addressed 
solely by considering the impedance and phase profiles. These include (i) the identi- 
fication of the mechanisms of amplification of the voltage response, (ii) the mecha- 
nisms of selection of the resonant and phase-resonant frequencies, (iii) the identifica- 
tion of the roles played by the intrinsic time scales and the time scales associated to 
the current inputs, (iv) the mechanisms that govern their interaction, (v) how all this is 
affected by changes in the model parameters, (vi) the relation between intrinsic STOs 
and subthreshold resonance, and (vii) the relation between subthreshold resonance 
and phase-resonance. 

Perhaps our intuition on the resonant properties of neuronal models is based on 
the dynamics of the so-called k-co systems (55)-(56) (in Appendix A.3), which were 
used to build the resonate-and-fire models introduced in [31]. These systems dis- 
play intrinsic oscillations with natural frequency £? nat = co for all values of A (58). 
They also exhibit resonance and phase-resonance. The natural, resonant and phase- 
resonant frequencies (59) coincide for A = 0 and they are different for other values 
of A. Although X-co systems have been used to investigate the dynamics of resonant 
neurons [31], they are not representative of linearized neuronal models [11]. They 
rather correspond to cases where the voltage and gating variables evolve with compa- 
rable rates (e = 1) (see (61)-(62) in Appendix A.3), and leave out the large class of 
neuron types that have a strong time- scale separation (e is small) in the subthreshold 
regime. 

The fact that resonance may occur in the absence of intrinsic STOs for non- 
negligible parameter regimes [11, 29] implies that resonance is not necessarily re- 
flecting the amplification of an existing intrinsic oscillation by an oscillatory input 
current. Instead, resonance is uncovering the ability of the neuron to operate at time 
scales that are neither intrinsic nor imposed by the oscillatory input, but they emerge 
as the result of the interaction between a neuron's intrinsic time scale (determined by 
the neuron's intrinsic properties) and the time scales of the oscillatory input. These 
emergent time scales are likely to be the ones that play a significant role in network 
interactions. 

To address these mechanistic issues from a geometric/dynamic perspective, we 
have extended the classical phase-plane analysis approach to include the effects of 
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oscillatory inputs. This approach consists of projecting the three-dimensional space 
(for v, w and t) onto the two-dimensional plane (for v and w) and viewing the pro- 
jection of the two-dimensional u-nullsurface (spanned by the u-nullcline for the au- 
tonomous system and t) as a moving one-dimensional u-nullcline as t progresses. 
Trajectories track the motion of the i;-nullcline and the fixed point with a speed that 
depends on the underlying vector field and the input frequency /. The shapes of 
the limit cycle trajectories depend on /. A system exhibits resonance if for some 
value of / the amplitude of the limit cycle trajectory in the v -direction is larger than 
for / = 0. Qualitative predictions on the effect of changes in parameters on both 
resonance and phase-resonance can be made by looking at the phase-plane and the 
associated envelope-plane diagrams. Approaches including "moving nullclines" have 
been used before to investigate the mechanisms of synchronization of neuronal [44- 
47] and other systems [48, 49] but they have not been used before in the context of 
the analysis of subthreshold and phase resonance. 

From a geometric perspective, we view both resonance and phase-resonance as 
the result of the interaction between the time scales of the neuron (determined by the 
geometry of the unperturbed phase plane and the intrinsic time- scale separation) and 
the oscillatory input (captured by the moving u-nullcline). The resonant frequency, if 
it exists, is the input frequency at which the voltage responds optimally to the oscil- 
latory driving current, which is captured by the cyclic movement of the f-nullcline. 
The trajectory's response is neither too fast nor too slow so that the voltage is able 
to reach a higher value than for other input frequencies. The time scale associated to 
the zero-phase frequency is also an emergent time scale, and is optimal in the sense 
that the trajectory is neither too fast nor too slow so that the trajectory intersects the 
moving f-nullcline at the exact time at which the latter reaches its maximum; i.e., 
both the voltage response and the input current peak at the same time. We showed 
that these two phenomena are captured by the shape of the envelope-plane diagrams 
not only for linear models, but also for nonlinear ones. 

The concept of time scale for neural models is not always easy to quantify. In some 
limiting cases the time constants provide reliable information about the rate of change 
of the participating variables. This is the case for the one-dimensional passive mem- 
brane equation and for the so-called slow-fast systems [47] . However, time constants 
do not always capture the effective time scales. For the purpose of our study, we have 
qualitatively characterized the effective time scales of the isolated neurons by looking 
at their response to instantaneous DC (tonic) inputs and the shapes of the correspond- 
ing trajectories. We have identified the effective time scales of the isolated neurons 
that are relevant for their interaction with oscillatory inputs with the time scales that 
govern the behavior of the initial portion of the trajectory from its initial point until 
it reaches its highest voltage value. Geometrically, this is determined by the crossing 
point with the u-nullcline displaced by A- m units (dotted-red v -nullclines). The effec- 
tive time scales defined in this way may or may not represent the dynamic behavior 
of the autonomous trajectories for all subsequent times t, depending on the parameter 
values, but they are appropriate to describe the interaction between neurons and oscil- 
latory inputs. As we showed, trajectories reverse direction after crossing the moving 
v -nullclines. 

This notion of effective time scales for the underlying autonomous system is con- 
nected with its eigenvalues. However, our discussion of resonance is independent of 
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the behavior of the autonomous system for large values of t . Therefore, the mecha- 
nisms of generation of resonance are in general independent from those responsible 
for the generation of intrinsic STOs. 

The mechanism of generation of intrinsic STOs in linear systems is well under- 
stood [50]. The natural frequency (J2 nat ) corresponds to the eigenvalues' imaginary 
part. The eigenvalues' real part affects the amplitude of solutions, and consequently 
the evolution of trajectories in the phase plane, without affecting the oscillation fre- 
quency. The initial time interval of transient behavior that determines the effective 
time scales depends on both the eigenvalues' real part and £? na t, and not only on £? n at- 
Only when the effect of the eigenvalues' real part is small enough, both the generation 
of intrinsic STOs and resonance are governed by a similar underlying mechanism as 
occurs for the k-co systems discussed above. 

To our knowledge, no previous study has addressed questions on subthreshold 
resonance in nonlinear systems from an analytical perspective. In this paper, we have 
extended our analysis for linear systems to include simple, piecewise-linear types of 
nonlinearities in both the equations for v and w. Our main question was: to what 
extent does the linear prediction capture the nonlinear effects? We found that when 
the i; -nullcline is nonlinear the differences between the nonlinear response and the 
linear prediction increase not only with increasing values of the input amplitude A- m 
but also with increasing levels of the time- scale separation between the voltage and 
the gating variable (decreasing values of e). These differences almost disappear when 
both equations evolve at comparable rates. In contrast, voltage responses are almost 
insensitive to nonlinearities located in the gating variable equation. In these two latter 
cases, the nonlinearities are present in the system but the voltage response does not 
detect them. More research is needed to understand whether these findings play out 
in more realistic nonlinear models. 

Resonance has been first studied in the damped harmonic oscillator subject to 
sinusoidal forcing (see Appendix A.4). Similar to the k-co system, the natural and 
resonant frequencies coincide when the system is undamped and they are different in 
the damped case. Unlike the neural models discussed in this paper, the over-damped 
system (real eigenvalues) does not exhibit resonance. The damped harmonic oscil- 
lator can be rewritten as a system of two first order equations (64)-(65). However, 
differently from the systems considered in this paper, the sinusoidal forcing is lo- 
cated in the second equation rather than the first. In other words, the forcing term 
does not directly affect the dynamics of the variable for which we compute the re- 
sponse to the oscillatory input (analogous to v), but its derivative (analogous to wj). 
Therefore, one should not necessarily expect resonance to occur in the absence of 
intrinsic oscillations. The geometric ideas developed in this paper can adapted to 
these systems. The moving nullcline will be the analog to the u; -nullcline (oblique 
with negative slope) instead of the v -nullcline. The latter will be horizontal and 
fixed. 
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Appendix A: Forced Two-dimensional Linear Systems: Eigenvalues, Natural 
Frequency, and Impedance and Phase Profiles 

We consider the following forced two-dimensional linear system: 

dx 

— = ax + by + A in sm(£?0> (41) 
dt 

dy 

— =cx +dy, (42) 
dt 

where a,b, c, and d are constant, Q > 0 and A- m > 0. 
A.l Intrinsic Oscillations and Natural Frequency 

The Jacobian matrix for the corresponding autonomous system {A- m = 0) is given by 

'-(:$)• <«> 

The roots of the characteristic polynomial are given by 

(a + d)± J (a - d) 2 + Abe 

7*1,2= ^ ' ^ ^ 

From (44), the unforced system displays oscillations for values of the parameters 
satisfying 

Abe + (a - df < 0, (45) 
with the natural frequency given by 

j-4be-(a-d) 2 

^nat = " " • (46) 

A. 2 Voltage Response to Sinusoidal Inputs: Impedance Amplitude and Phase 
System (41)-(42) can be written as a second order linear ODE, 

d 2 x dx 

— T -(a + d) Y{ad- bc)x = A m \-d sin(X2f) + Q cos(^0l- (47) 

dt 1 dt L J 

The particular solution to system (47) has the form 

x p (t) = A ou t,i sin(^0 + A ou t,2cos(^0 (48) 
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with 



and 



(ad -be- Q 2 )d + (a + d)Q 2 A 
AouU " [ad-bc-Q2? + (a + d)2a* Ahl m 



_ (ad -be- Q 2 )Q - (a + d)Qd 
AouU = [ad-be-C2 2 ] 2 + (a + d) 2 (2 2 A ' m ' m 
The solution (48) can be written as 

x p (t) = A 0Ut sin(Qt -0). (51) 

The impedance amplitude and phase are given by (see Appendix A.5) 

7 2 (G\ — - A out,i + A out,2 d 2 + Q 2 

K } ' A? A? M - be - Q 2 ] 2 + (a + d) 2 Q 2 K } 

and 

_i/A 0Ut , 2 \ _! (ad-*c-tf 2 )tf-(a + </)tf</ 

= — tan = tan — , (53) 

^ VAout.1/ (ad-bc- Q 2 )d + (a + d)Q 2 

respectively. The impedance Z peaks at the resonant frequency 



^ res = y - J2 + yz?2 c 2 2 ^cd - 2d 2 Z?c, (54) 
provided the quantity inside the radical is positive. 
A. 3 The X-co Systems with Sinusoidal Forcing 

The so-called X-co systems [51] with sinusoidal forcing have the form 

dx 

— = — Xx — coy + B[ n sin(£?0> (55) 
dt 

— = cox — Xy, (56) 
dt 

with X > 0, co > 0 and #i n > 0. System (55)-(56) can be written as the following 
second order linear ODE: 

^-4 + 2A— + (X 2 + <w 2 )jt = £iJA sin(X?0 + cos(Qt)]. (57) 
dt 1 dt v 7 L J 

The eigenvalues and natural frequency of the autonomous system are given by (see 
Appendix A.l) 



n 2 = — X ± \l -co 2 and £? nat = (58) 
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The resonant and phase-resonant frequencies are given by (see Appendix A. 2) 



^ res = y~A 2 + co^X 2 + co 2 and Q phdS = ^ co 2 - X 2 , (59) 
provided the quantities inside the radicals are positive. For A = 0, 

^nat = ^res = ^phas- (60) 

System (55)-(56) can be transformed into a system of the form (14)— (15) by defin- 
ing 

v = Xx, w = coy, i = Xt, (61) 

and 

co 2 

A.4 The Harmonic Oscillator with Sinusoidal Forcing 

The equation for the harmonic oscillator with sinusoidal forcing reads 

d 2 x dx 

~dt 2+P li + yX = ° m sin( ^ } ' (63) 

where ft > 0 and y > 0. By defining y = —dx/dt, (63) can be rewritten as the fol- 
lowing two-dimensional linear system of ODEs: 

dy 

= yx - fly - C in sm(Gt). (65) 

dt 

Differently from the systems considered in this paper, and the general form (41)- 
(42), the sinusoidal forcing is located in the second equation rather than the first. 
We note that the two equations are not interchangeable since we are following the 
convention that the first equation describes the dynamics of the variable (x) for which 
we compute the response to the oscillatory input. 

The eigenvalues of the autonomous system are given by 

n,2 = - • (66) 

The natural frequency of the autonomous system is given by 



^nat= " = \ Y ~~^' (67) 
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provided 4y — /3 2 > 0. The impedance is given by 

1 

Z(Q) = ^ =-^r. (68) 

v } (y -n 2 ) 2 + p 2 n 2 v } 

Note that the formula (54) is not applicable in this case and the impedance has to be 
calculated separately. The impedance peaks at 



^res = yK-y, (69) 

provided the quantity inside the radical is positive. For /3 = 0, £? nat = £? res = y/y. For 
other values of /3 for which both £? nat and £? res are defined, £? nat ^ £? res . If the eigen- 
values are real, the over-damped harmonic oscillator does not exhibit resonance since 
the inequalities Ay — fi 2 < 0 and 2y — fi 2 > 0 cannot be simultaneously satisfied. 
System (64)-(65) can be rescaled by defining 

t = /3t and y = -y (70) 

y 

and substituting into system (64)-(65). The resulting equations read 

dx y 



dt 



2 y, (71) 



^ = x - y - — sm(t2i/P). (72) 
dt y 

A. 5 Oscillatory Inputs: Additional Calculations 

For a sinusoidal input of the form F(t) = A- m sin (£2t) the system's output will be a 
function 

X(t) = A out sin(Qt-(t>). (73) 
Equation (73) can be rewritten as follows: 

X(t) = A out cos(j)sm(f2t) — A ou t sin 0 cos (£2t) (74) 



or 



with 



X(t) = Aout.1 sm(Qt) + A out ,2cos(^0 (75) 



Aout,i = Aout cos 0, A ou t,2 = - Aout sin 0. (76) 
Solving for A ou t and 0 we obtain 



^out — ^out,l + ^out,2 (77) 
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and 



0 = — tan 



-l 




(78) 



From (77) 



Z 2 {Q) = 




out, 2 



(79) 



Appendix B: Biophysical Ih + /Nap and /ks + /Nap Models 

B.6 4 + 7 Nap Model 

This model has been proposed in [37]. It has a persistent sodium current and a two- 
component (fast and slow) h-current given by 7 p = G v p(V — E^a) = Gp/?oo(^0 x 
(V — #Na) and 7 n = G^riV — £ n ) = G n (c/r/- + c 5 r 5 )(y — £h), respectively, with 
£"h = —20 mV, #Na = 55 mV, Cf = 0.65, and c 5 = 0.35. The voltage-dependent ac- 
tivation/inactivation and time constants are given by 



p^(V) = l/(l+e-™^), 
t p (V) = 0.15, 

r / , 00 (y ) = i/(i + ^+ 79 - 2 )/ 9 - 78 ), 

T r/ (V) = 0.51/( e ^- 1 - 7 )/ 10 + e -^ +340 >/ 52 ) + 1, 

r s , 00 = l/(l+^+ 71 - 3 )/ 7 - 9 ), 
Trs (V) = 5.6/{e {V - lJ ^ u + e "^+ 2 «»/43) + 1. 



For the two-dimensional model used in this paper, the Cf = 1 and c s = 0. 
B.7 7 Ks + /Nap Model 

This model has been proposed in [37]. It has a persistent sodium current and a slow 
potassium (M-type) current given by / p = G v p(V — E^ a ) = G v poc(V)(V — E^ a ) 
and /ks = G q q(V — E Xk ) with E^a = 55 mV and E^ = — 90 mV. The voltage- 
dependent activation/inactivation and time constants are given by 



^(7) = 0.15, 

4oo(V) = l/(l+e 

^r(V) = 90. 



(V+38)/6.5 



(V+35)/6.5 
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